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Abstract 

In this paper we discuss a closed-form approximation of the likelihood functions of 
an arbitrary diffusion process. The approximation is based on an exponential ansatz 
of the transition probability for a finite time step At, and a series expansion of the 
deviation of its logarithm from that of a Gaussian distribution. Through this procedure, 
dubbed exponent expansion, the transition probability is obtained as a power series in 
At. This becomes asymptotically exact if an increasing number of terms is included, 
and provides remarkably accurate results even when truncated to the first few (say 3) 
terms. The coefficients of such expansion can be determined straightforwardly through 
a recursion, and involve simple one-dimensional integrals. 

We present several examples of financial interest, and we compare our results with 
the state-of-the-art approximation of discretely sampled diffusions [A'it-Sahalia, Journal 
of Finance 54, 1361 (1999)]. We find that the exponent expansion provides a simi- 
lar accuracy in most of the cases, but a better behavior in the low-volatility regime. 
Furthermore the implementation of the present approach turns out to be simpler. 

Within the functional integration framework the exponent expansion allows one to 
obtain remarkably good approximations of the pricing kernels of financial derivatives. 
This is illustrated with the application to simple path-dependent interest rate deriva- 
tives. Finally we discuss how these results can also be used to increase the efficiency 
of numerical (both deterministic and stochastic) approaches to derivative pricing. 



1 INTRODUCTION 

Continuous-time diffusion processes is the basis of much of the modeling work performed 
every day in Finance and Economics, from portfolio optimization and econometric appli- 
cations, to contingent claim pricing. Indeed, since Bachelier's 1900 doctoral thesis on the 
dynamics of stock prices [I], many economic variables subject to unpredictable fluctuations, 
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have been modeled by stochastic differential equations of the form 

dY t = v v {Yt)dt + a y (Y t )dW t . (1) 

Here n y (y) is the drift, describing a deterministic trend, and o~ v (y) > is the volatility 
function, describing the level of randomness introduced by the Wiener process (i.e., white 
noise) , dWt . The main reason for the popularity of this class of models is probably that in 
continuous time one can perform analytic calculations using the instruments of stochastic 
calculus, and the powerful framework of partial differential equations. 

In particular, for the few cases for which the process fT]) is exactly solvable, one can 
derive closed-form solutions for the associated transition probability. The latter contains 
all the statistical properties of the financial quantity modeled by the diffusion, and can be 
exploited in a variety of ways, including the derivation of no-arbitrage prices for financial 
derivatives in complete markets. The milestone results derived by Black, Scholes and Merton 
[2j [3] , Cox, Ingersoll and Ross [4 , or Vasicek [6] , are among the most significant examples of 
the amount of progress in Economics that has been done using integrable continuous-time 
diffusion processes. 

Nonetheless, an accurate description of the market observables requires in general more 
sophisticated models than those for which an analytic solution is available. These are usually 
tackled by means of numerical schemes ultimately relying on a discretization of the diffusion, 
obtained by replacing the infinitesimal time dt with a finite time step, At. These approaches 
typically involve either solving numerically a Kolmogorov partial differential equation, or 
a Monte Carlo sampling of the diffusive paths. The approximate results obtained in this 
way become exact only approaching the limit At — > 0, and this can be done with some 
computational effort. In addition, a drawback of these numerical methods, more specific to 
econometric applications, is that they do not produce closed-form expressions for the tran- 
sition probability. These are crucial for maximum-likelihood estimations of the parameters, 
say 9, of model diffusions. In fact, economic data are generally available on discrete sets of 
observations, usually well spaced in time, say weekly or monthly. As a result, only if the 
transition probability, p(Xn + x\^ t , At\Xi\ t ] &), associated with each time interval At of the 
series t — I At (I = 1, . . . n), is known in closed form, the maximum likelihood function, 

1 " 

= - E lo SP( X (l+i)At, IMXiau 9) , (2) 
i=i 

can be analytically maximized over 9. If this is not the case, one has to repeat the numerical 
calculation of the transition probability for every value of 9 needed for the determination of 
the maximum, e.g. by means of an optimization algorithm. This can be clearly very time 
consuming. 

Motivated by this difficulty, A'it-Sahalia recently proposed a method to approximate the 
transition probability for one dimensional diffusions [7] by means of a Hermite polynomial 
expansion. This was applied to a variety of test cases, and the accuracy of the method was 
clearly demonstrated. 

In this paper, we utilize the exponent expansion - a technique introduced in chemical 
physics by Makri and Miller [5] - to derive a closed-form short-time approximation of the 
transition probability of the diffusion process (|TJ) . The aim is to obtain an analytic approxi- 
mation which is as accurate as possible for a time step At as large as possible. On one hand, 
this allows one to derive approximations of financial quantities that are very accurate even 
for sizable values of the time step, and to derive closed-form expression for the maximum 
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likelihood function @. On the other, it allows a reduction of the computational burden of 
numerical schemes as the limit At — > can be achieved with larger time steps, i.e., with 
less calculations. Our approach is similar in spirit to the one of Ref. [7J, reviewed in Section 
12.11 but it overcomes some of its shortcomings, and it is of simpler implementation. In 
particular, the coefficients of the exponent expansion can be expressed in terms of one di- 
mensional integrals that can be easily calculated numerically. In addition, we show how we 
can apply our approach to any sufficiently regular volatility function, even when the latter 
does not have an analytic expression but it is specified numerically through an interpolation 
procedure, as it is the typically the case for local volatility models. 

The possibility to use Makri and Miller's technique to derive approximations of the 
transition probability was originally hinted by Bennati et al. in Ref. [HI HO]- Here we 
explore this route, giving derivations for a generic diffusion process with state-dependent 
drift and volatility, and we study the reliability of the exponent expansion by applying it to 
several test problems of financial interest. 

Through the exponent expansion, the transition probability is obtained as a power series 
in At which becomes asymptotically exact if an increasing number of terms is included, and 
provides remarkably accurate results even when truncated to the first few (say, n = 3) 
terms. Two derivations are offered, the first by means of Kolmogorov's forward equation 
[TT] (Sec. 12. 2p . and the second introducing a slightly different formalism (Sec. I2.2.1[) . The 
latter, once the problem is formulated in terms of Feynman's path integrals [IHEI], allows 
the generalization of the exponent expansion to the calculation of the pricing kernel of 
financial derivatives whose underlying follows the considered diffusion. This allows in turn 
the derivation of simple approximations for the price of such contingent claims (Sec. [3]). In 
Sections 12.31 and 13. 3[ we illustrate the exponent expansion through the application to the 
Vasicek, the Cox-Ingersoll-Ross, and the Constant Elasticity of Variance models, and in 
Section [4] we discuss its application to Monte Carlo and deterministic numerical methods 
within the path integral framework [TH UJ3 HH H7J [TSl [TU] . Finally, we draw our conclusions, 
and we discuss future developments in Section [6l 



2 TRANSITION PROBABILITIES OF DIFFUSION 
PROCESSES 

Let us consider the problem of estimating the transition probability, p y (y, At\yo), associated 
with one-dimensional continuous-time diffusion processes of the form (|TJ) . This represents 
the likelihood that the random walker following the process (UJ ends up in the position y at 
time t = At, given that it was in y at time t = 0, and satisfies the Kolmogorov forward (or 
Fokker-Planck) equation [TT] : 



d t py(y,At\y ) 



-d v ny{x) 



d 2 y <j y {y) 2 



py{y,Myo) 



(3) 



In this Section we will consider two short-time approximations of the transition probability 
above that can be derived in closed form. We will start by reviewing the Hermite poly- 
nomial expansion, recently introduced by Ait-Sahalia [7], and then describe the exponent 
expansion. 
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2.1 Review of the Hermite polynomials expansion 

The first step of A'it-Sahalia's derivation is to transform the original process in an auxiliary 
one, say X t , with constant volatility a x . Following Ref. [7], this can be achieved in general 
through the following integral transformations 

X t = 1 (Y t )^±o x T (4) 
J cr v {z) 

where the choice of the sign is just a matter of convenience depending on the specific problem 
considered. The latter relation defines a one to one mapping between the X t and Y t processes 
as the condition cr y (z) > ensures that the function x = 7(2/) defined by (|4]) is monotonic, 
and therefore invertible. Q A straightforward application of Ito's Lemma [IT] allows one to 
write the diffusion process followed by Xt as 

dX t = fjL x (X t )dt + <r x dW t , (5) 

with 



Hx{x) = ±a x 



Hx)) _ I da y ! . 
a y ( 7 -i(x)) 2 dy 17 {X) \ 



(6) 



where y = ~f~ 1 (x) is the inverse of the transformation (H|). 

Using Hermite polynomials, it is possible to show [7 , that a short time approximation 
of the transition probability can be expressed up to order N as 

dzfj, x (z))^2c n (x,x ) — , (7) 

-'0 n=0 

where <j>{x) is a standard normal distribution. Here the coefficients c n (x,xo) can be derived 
by solving the recursive equation 

c n (x,x ) = n(x - x y n [ dz(z - x ) n ~ l ( X(z)c n ^ 1 (z 7 x ) + \d 2 z c n ^i{z , x ) J , (8) 



Jx 

for n > 0, with cq(x,xq) = 1, and 
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X(x) = - i (^(a;) 2 + QrMs)) • (9) 

Finally, the transition probability for the process Y t can be determined through the 
Jacobian of the transformation ([4]) giving 

/ A , 1 \ Px(j(y), At\x ) 
Py{y,At\yo) = <?x r-s • (io) 

The expansion above turns out to provide very accurate results in a variety of test cases 
for which the exact expression of the transition probability is available. Some of these 
examples will be considered in detail in the following when we will compare the results 
obtained with the approximation Eqs. (|7|)- (|10p . with those given by the exponent expansion 
described in this paper, and introduced in the following Section. 



^^For a discussion of the regularity conditions on the drift and volatility functions see e.g., Ref. [7]. 
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2.2 The Exponent Expansion 

In this Section, we derive the exponent expansion for the transition probability for the 
process ([T]), p y (y, At\y ). In order to make the derivation easier, it is convenient to transform 
the original process in the constant volatility one Eq. ([5]) by means of the transformation 
©■ 

For At — > the transition probability for X t is dominated by the diffusive Gaussian 
component and therefore reads: 



p x (x, AT\x ) 



1 



2wa 2 r At 



1/2 



exp 



(x - x f 
2a? At 



(11) 



Guided by this observation, in order to find an expression for the transition probability 
associated with Eq. ([5]) which is accurate for a time At as long as possible, we make the 
following ansatz: 



p x (x,At\x ) 



1 



: exp 



(x - x Q ) 2 
2(7? At 



W(x,x ,At) 



y^alAt 

Such transition probability must satisfy the Kolmogorov forward equation [11] : 

f 



d t p(x, At\x ) 



-d x p x {x) + ^<7 2 x d 2 x 



p x (x, At\xo) 



(12) 



(13) 



Note that since the auxiliary process X t has a constant volatility, the latter equation is 
mathematically much simpler than the one for the process Y t ([3]). Equation (|13p implies in 
turn that the function W(x, xq, t) satisfies the relation: 



8 t W = -p x d x W + \a 2 x d 2 x W - X -al (d x W) 2 + d x p x - 
Expanding W(x, xo, t) in powers of At, 

oo 

W(x,x ,At) = Y, W ^ x > x °) Atn ' 



d x W+^ 



(14) 



(15) 



substituting it in Eq. (|14p , and equating equal powers of At leads in a straightforward way 
to a decoupled equation for the order zero in At giving 



1 f x 

W (x,x ) = j / dzp x (z) , 

® x J xo 

and to the following set of recursive differential equations: 

1 



(n + 1)W, 



n+l 



-[x - x )d x W n+1 + 



^ m—n 

n a l X! 9 xW m d x W n ^ m + 8 nfi d x p, x 



m=Q 



(16) 



(17) 
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In particular, for n = 0, 1, 2 Eqs. (jTTJ) read: 

W 1 (x,x ) = -(x - x )d x W 1 (x,x ) + 



-^-p\i x (x) 2 + ~d x fj, x (x) 



2W 2 {x, x ) 
3W 3 {x, x ) 



-(x - x )d x W 2 (x,x ) + ^aldlW 1 (x,x ) , 
-(x - x Q )d x W 3 (x,x Q ) + ^a 2 x d 2 x W 2 (x,x Q ) , 
-aKdxWt&xo)) 2 . 



(18) 
(19) 

(20) 



The differential equations above (|17[) are all first order, linear and inhomogeneous of the 
form 

nW n (x, x ) = -(x - x )d x W n (x,x ) + A n -i(x,x ) , (21) 

where A n _i(x, Xq) is a function that is completely determined by the first n — 1 relations. 
It can be readily verified by substitution and integration by parts that the solution of (|2ip 
reads 



W n (x,xo) 



^r _1 A„-i(x + (x - xo%x ) 



This, for n — 1,2,3, after some manipulations, gives: 



Wi{x,x ) 
W 2 (x,x ) 
W 3 (x,x ) 



1 

Ax 



dzV e s(z) 



2Ax 2 



^ [V eS {x) + V eS (x ) - 2^1(1,10)] 



2Ax 4 
3a? 



Ax / dzV eS {z) 2 



dzV s(z) 



^W 2 (x,x ) + ^3 [d x V cS (x) - d x V cS (x )} 



(22) 

(23) 
(24) 

(25) 



where Ax = x — xq , and, for reasons that will be clearer in the next Section, we have also 
introduced the 'effective potential as the following quantity with dimension time -1 : 



V cS (x) = —^^ x {x) 2 + \d x ^ x (x) 
2ai 2 



(26) 



Note that, from Eq. ([9]), V c s (x) = — 2X(x). The first order correction can be rewritten as 



1 f At 

W 1 (x,x ) = —J dtV cS (x Q + tAx/At) 



(27) 



leading to the interpretation of this term as a time-average of the effective potential over 
the straight line, constant velocity (Ax/ At) trajectory between Xq and x. Similarly, the 
leading term in W 3 {x, Xq) (i.e., the one proportional to the lowest power of the volatility) is 
proportional to the variance of the effective potential over the same trajectory. Note that 
the corrections W n (x, Xq) are well defined in the limit Ax — > 0. In particular, for n — 1,2,3 



(i 



it is not difficult to show that 



lim Wx{x,x Q ) = V eS (x ) , (28) 

2 

lim W 2 (x,x ) = ^d 2 x V cS (x), (29) 

lunW 3 (x,x ) = -%L( dx V eS (x)f + ^f-d$V eS (x). (30) 



Finally, the transition probability of the original diffusion |T]) is recovered by means of 
the Jacobian transformation ([TU]) . 

The form of the trial transition probability represents the main difference of the present 
approach to the one described in Section T2.1 1 Ref. [7J, which is otherwise very similar in spirit. 
In fact, the latter expands in powers of At the exponential exp [— W(x, xq, At)] rather then 
just the exponent, as we do here, instead. As it will be shown explicitly in the following, the 
present choice gives rise to a distinct approximation scheme for n > providing generally 
a similar level of accuracy but remarkably simpler mathematical expressions. In particular, 
all the calculations related to the test cases presented in the following Sections have been 
easily obtained by hand, i.e., without the aid of symbolic calculation computer programs. 
This is because, by keeping the exponential form of the ansatz, one formulates a guess 
which is closer to the exact one. The latter, can be expected to have an exponential form in 
order to satisfy the Chapman-Kolmogorov property of Markov processes [H] . In addition, 
the exponential choice of the ansatz automatically enforces the positive definiteness of the 
transition density which is not granted in the approach of Ref. [7J. In fact, as it will be 
shown explicitly in Sec. 12.31 in the limit of very small volatility (a x — > 0), when the effect 
of the noise disappears, and the transition density converges to a Dirac's S distribution, 
the expansion of Ref. [7J breaks down as the transition probability becomes negative. On 
the contrary, the exponent expansion remains well defined and accurate also in this limit. 
Indeed, the first terms of the expansion in At can be also derived through a small volatility 
expansion of the transition density, as it will be discussed in Sec. 13.11 

Similarly to the approximation developed in Ref. [7J, the exponent expansion has in 
general a finite convergence radius which is a decreasing function of the volatility. As it 
will be shown in the following, for the values of volatilities and At relevant for financial 
applications the exponent expansion turns out to be very accurate even when truncated to 
the first few terms. 

2.2.1 Alternative derivation 

The term Wq{x, Xq) in the exponent expansion is somewhat different from the higher order 
terms. In fact, it is defined by Eq. p6[) which is decoupled from the recursive system (jTTJ) . 
Indeed, it is possible to obtain the same result for the exponent expansion by expressing the 
transition density as 

p x (x,At\x ) = e- w °( x '*°^ Px (x,At\x ) , (31) 

and looking for an approximate expansion of the form (|12j) for $ Pa ,(a;, At|xo). It is easy to 
show by direct substitution in the forward Kolmogorov equation (|13| that & Px (x, At\xo) is 
the solution of 

H x $ Px (x, At\x ) = -d t $ Px {x, At\x ) , (32) 
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where H x is the "Hamiltonian" differential operator 



-d 2 x + V cS (x) 



(33) 



and V e g(x) is the effective potential of Eq. ([26)) . As a result one can equivalently derive the 
exponent expansion by substituting in Eq. (|32|) the following trial function 



*p« (x,At\x ) 



exp 



2a 2 . At 



(34) 



which does not contain the term Wo(x,xq). This observation will be used in Section 3 to 
generalize the exponent expansion to the pricing kernel of financial derivatives. 



2.2.2 Exponent Expansion Coefficients in terms of the the Original Diffusion 

The derivations of the exponent expansion, and of the Hermite polynomial approximation 
of Section [2TTT both rely on the introduction of the auxiliary process ([5]), and the integral 
transformation ([4}. As a result, the expressions obtained are easy to handle if the the 
volatility of the original process a y (y) is such that both the function 7(2/) (j4|) and its inverse 
7 _1 (a;) admit a closed-form expression. In fact, in this case the effective potential (f26|) 
- or the function A(z) - has a closed-form expression, and the determination of the 
coefficients of the expansion can be determined either analytically or through numerical 
quadrature. This is the case for the examples considered in the following. However, these 
are very special cases. In fact, very few volatility functions have a reciprocal for which 
a primitive exist. On the contrary, in many practical applications, e.g. in local volatility 
models jS] , the volatility is specified numerically through a fit to an volatility function. For 
all these cases the application of the exponent expansion or the approach of Section 12.11 
becomes cumbersome and computationally demanding, because it requires the numerical 
inversion of the integral equation (j4|) . 

However, at a more careful analysis, it turns out that it is possible to circumvent this 
difficulty, and to eliminate any dependence on the function 7 _1 (cc) in the expressions for the 
exponent expansion. This makes its application straightforward for any diffusion process, 
irrespective of the analytic tractability of the specified volatility function. The first step to 
do this, is to note that, according to the Jacobian transformation (p~0|) . x and xo in Eq. (fP2|) 
need to be calculated for 7(2/), and 7(2/0): respectively. As a result, one can express the 
exponent expansion as 



Py(y,&t\ya) = 



with the notation 



1 



: exp 



(i(y) - i(yo)) 2 

2At 



-J2Wn(y,y )Af> 



n=0 



G(y,yo) = G(x,a;o)U=7(y),x=7( 1 ,o)) 



(35) 



(36) 



and taking, without any loss of generality a x = 1. Now, using Eq. (|16p . by means of the 
change of variables y — r y~ 1 (x) one can express the zeroth order term as 

f v dz 

W (y,y ) = - — — fi x {z) 



yo u v 



(Ty(z) 
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where, using Eq. 



(Jy(z) 2 dy 



(z) 



(37) 



As anticipated, the expression for Wo does not contain any reference to the function j^ 1 (x) 
so that it can be calculated even if there is no closed-form expression for 7(2/) = dz/<r y {z) 
available, or if the latter is not analytically invertible. Similarly, using Eqs. ([23" |) -(f25 f one 
can easily find 



Wi{y,y ) 
W 2 (y,y ) 
W 3 (y,y ) 



1 

A£ 
1 



dz 



2Ay 2 
1 



V cS (y) + V c s(yo)-2W 1 (y,y ) 



2Ay 



f,4 



Ay 



-^W 2 {y 7 y 



dz 

VO a v( Z ) 



V eS (zf 



dz 
yo a y( z ) 

dyVes(y) - d y V e g(y) 



V eS {z) 



Ay 2 "W"' 4Ay 3 
where Ay = j(y) — 7(2/0)1 and the effective potential now reads 



(38) 
(39) 

(40) 
(41) 



As a result, the exponent expansion can be easily applied for any specification of the volatility 
function <J y (y). Indeed, given any drift and volatility function one can immediately calculate 
- possibly numerically - their derivatives and the function 7(2/). From these quantities, 
one then easily calculates the effective potential (|4"T1) . and the coefficients of the exponent 
expansion ([55 ]) - ([ID"] ) . 



2.3 Examples 

The application of the exponent expansion to a generic diffusion process of the form (fT|) 
is rather straightforward and reduces to the calculation of one dimensional integrals. In 
this Section, we illustrate this procedure for a few test cases, namely for the Vasicek [6], 
the Cox-Ingersoll-Ross [3], and the Constant Elasticity of Variance [Hj diffusion processes. 
We will compare the results of the exponent expansion with the exact results available in 
literature, and with the approach of Ref.[7]. 



2.3.1 Vasicek diffusion 



We first consider the Ornstein-Uhlenbeck diffusion proposed by Vasicek [5] as a model for 
the short-term interest rate: 



dX t = a{b - X t )dt + adW t 



(42) 



where a, b, and a are positive constants representing the mean-reversion level, the veloc- 
ity to mean reversion, and the volatility, respectively. This model is integrable and the 
corresponding probability density function is Gaussian: 



p ex (x, At\x ) 



1 

(2 7 ra 2 ) 1 /2 



exp 



[(x 



-aAt 



2(7 2 



a)} 



2 1 



(43) 
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Figure 1: Accuracy of the exponent expansion for the transition density of the Vasicek model 
(|42f . for a = 0.0717, b — 0.261 and xq — 0.1. Panel (a): comparison between the exponent 
expansion (stars) and the approximation of Ref. [7] (squares) for a = 0.02237 and At = 0.5. 
At order zero the two schemes are identical. The uniform error of the Euler approximation 
is also reported for comparison. Panel (b): comparison between the exponent expansion 
(stars) and the approximation of Ref. 7j (squares) in the regime of low volatility (<r = 0.01), 
for At = 0.5. 
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Figure 2: Accuracy of the exponent expansion for the transition density of the Vasicek model 
(|42]>, for a = 0.0717, b = 0.261 and x a =0.1. Maximum relative error for a = 0.02237 as a 
function of At: n = (continuous), n — 1 (dotted), n = 2 (long dashed) and n = 3 (short 
dashed). The inset is an enlargement of the 5-10 years region. 
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with 



a = a 



-2aAt 



2a 



(44) 



The exponent expansion of the Vasicek model can be easily derived using Eqs. (|16p . and 
(|23ti25|) with the effective potential, Eq. (|26j). 



and gives, 



p x (x, At\x ) 



V eS {x) 



1 



a 2 {b-xf 



2a 2 



a 
2 



(45) 



: exp 



(x - x ) 2 



W (x,x ) 



V2TT<r 2 At L 2a 2 At 
- W 1 (x,x )At-W 2 (x,x )At 2 - W 3 (x ) x )At 3 ] 

up to the third order in At (n = 3). Here 

a(x — b) 2 — a(xQ — b) 2 



Wo(x,xo) 
Wi(x,x ) 
W 2 (x,x ) 
W 3 {x,x ) 



GfT' 



2a 2 

:((x - b) 2 + (x - b) 2 + (x- b)(x - b)) - 
i 

— [V cS (x) + V cfi (x Q )-2W 1 (x,x )] , 



(46) 

(47) 
(48) 
(49) 



2Ax 3 



20ct- 



■i(x~b) 5 -(x -by- 



^Ax--^[(x-b) 3 -(x -b) 3 } 



2A 



6a 

^-(W^x^ojf 



3o 



A. 



^W 2 {x,x ) 



4Aa 



(50) 



with Aa; = x — xo- ft is interesting to note that the approximate transition probability 
obtained with the present approach reproduces exactly the expansion of the exact transition 
density Eq. (|43[) at the same order. 

On the other hand, the first two coefficients of the Hermite polynomials expansion (JT]) 
as quoted in Ref. [7] read: 



Ci (x,x ) 



c 2 {x,x ) 



1 



a(3b a — 3(x + xo)baa 

(—3 + x 2 a + xxoa + x^a)a 2 ) 
1 

36^ 



(51) 



a 2 (%V - 18xb 3 a 2 a + 3b 2 a{-6 + 5x 2 a)a 2 

- 6xba(-3 + x 2 a)a 3 + (3 - 6x 2 a + x 4 a 2 )a 4 

+ 2aer(— 36 + xa)(3b 2 a — 3xbaa + (—3 + x 2 a)a 2 )xo 

+ 3acr 2 (5b 2 a - Axbacr + (-2 + x 2 a)a 2 )xl + 2a 2 a 3 (-3b + xa)x 3 



(52) 



The fast convergence of the approximation scheme is illustrated in Figs. [J^l Here the 
percentage error of the exponent expansion with respect to the exact result f4"3"|) is plotted 
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for various At, and compared with the approach of Ref. [7]. The parameter choice, also 
taken from Ref. [7] corresponds to a sensible parameterization for interest rate markets. We 
adopt one year as unit of time, and we express the various parameters in this unit. The 
Euler approximation 

i a+i ^ / 1 ( ( x ~ x ° ~ Vx{xa)At) 2 \ 
p E (*,M\x ) = ^/^^exp { ^ j , (53) 

is also reported for comparison. The inclusion of each successive order allows one to increase 
dramatically the accuracy of the approximation so that the third order expansion has basi- 
cally a negligible error even for a sizable time step of order 6 months. Remarkably, for the 
considered example, the third order expansion allows one to estimate the 10 years transition 
probability with a relative error of less than 10 basis points (Fig. [2]). As illustrated in Fig.[l|- 
(a), in the present case the approach of Ref. [7] provides a slightly poorer level of accuracy 
for n < 2. In addition, it generally produces more complicated mathematical expressions. 
Furthermore, as shown in Fig.[T]-(b), in the regime of small volatility the exponent expansion 
still provides accurate results while the performance of the approach of Ref. [7] degrades. 
In fact, as anticipated, in the limit of small volatility (a x < 0.5%) the first order correction 
of the latter approach produces a negative transition probability signaling a break down of 
the scheme. 



2.3.2 Cox, Ingersoll and Ross diffusion 

The Vasicek model is probably too easy of a test case as the associated transition probability 
is Gaussian. In fact, since the exponent expansion has a leading term which is Gaussian, the 
higher powers in At just have to renormalize its average and variance in order to reproduce 
the exact result. It is interesting therefore to test the accuracy of the exponent expansion 
for a diffusion process that, while still integrable, generates a non-normal transition density. 
This is the case for the Feller's square root process [20] 

dY t =a{b-Y t ) + a y ^Y t dW t (54) 

which is the basis of Cox-Ingersoll-Ross model for the instantaneous interest rate [I]. The 
exact transition probability is given by [3] 

p ex (y, At\y ) = ce-("+^ (~) § I q (2^) , (55) 

where c = 2a/[ijy{l — exp (— a At))], q = 2ab/a 2 — 1 > 0, u = cyo cxp (— aAt), v = cy and I q 
is the modified Bessel function of the first kind of order q [5T] . 

As explained in Sec l2.2l since the volatility is not uniform, it is convenient to introduce 
the auxiliary process defined by Eq. (UJ), as X t = 7(Jt) = 2\fYtja y . The X t process follows 
Eq. with cr x = 1 and 

(i x {x) = -x , (56) 

x 2 

and q = q + 1/2. In this case the effective potential reads: 

V eS (x) = ~» x (x) 2 - -^- a - , (57) 
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and the first four terms of the exponent expansion in Eq. (|15p are: 

x a 



Wq(x,x ) 
Wi(x,x ) 

W 2 (x,x ) 
W 3 (x,x ) 



-q\og— + -(x 2 -xf i ) 
x Q 4 



1 



2Aa 



12 v 
1 



fj-x(x) - fJ- x (xo) - q 2 



1 1 

X Xo 



2Ax 2 
1 



2Ax 2 
3 



-x ) - aq [x - X ) 
[V cB (x) + V cS (x ) - 2W 1 (x,x )] 
-(W 1 (x,x )) 2 

[dxV e s(x) - d x Ves(x )} 



G(x) - G(x ) 



Ax 



W 2 (x,x ) + 



1 



Ax 2 K 1 " 4Ax 3 
where Ax = x — xq, d x V c s(z) — q~(l — q)/z 3 + a 2 z/A, and 



G(z) 



25 ^ P , 2 , r, a 

-a z - - — r + 7 z + 2apz - 
5 3 z A 



-ajz 



(58) 

(59) 
(60) 

(61) 
(62) 



and a = a 2 /8, (3 = q(q — l)/2, 7 = —a(q + l)/2. Finally, going back to the original process 
with Eq. (fit)]) , the transition probability reads: 



Py(y, At|y ) = p x (2 v fy/a y , At\ 2y/%/a y )/a yy /y 



(63) 



The coefficients of the expansion in Hcrmite polynomials Eq. as quoted in Ref. 7 
read instead: 



ci(x,x ) 



c 2 {x 1 x ) 



1 



m 2 a 2 - 486acr 2 + 9ct 4 



24xxot 4 

2 2/ 0^7 1 2 2\ 1 2242, 243 

xa a [—24:0 + x a )xq + x a a Xq + xa a Xq 

' r 9(2566 4 a 4 - 5126 3 a 3 cr 2 + 2246 2 a 2 o- 4 + 32bacr 6 - 15er 8 ) 



(64) 



576x 2 x£ 

+ Qxa 2 o 2 {-2Ab + a; 2 cr 2 )(16bV - 166acr 2 + 3a 4 )x 

+ x 2 a 2 cr 4 (672b 2 a 2 - 486a(2 + x 2 a)a 2 + (-6 + x A a 2 )a A )x 2 

+ 2xa 2 a A (A8b 2 a 2 - 24ba{2 + x 2 a)a 2 + (9 + x A a 2 )cr A )xl 

+ 3x 2 a 4 CT 6 (-16fe + x 2 a 2 )x% + 2x 3 a 4 a 8 xl + x 2 a 4 a 8 x% 



(65) 



The accuracy of the exponent expansion in this case is illustrated in Figs. [3] and [U 
Similarly to the case of the Vasicek diffusion, the exponent expansion is characterized by a 
remarkably fast convergence by including successive terms of the approximation so that n = 
3 provides already a virtually exact representation of the transition density, for At ~ 1 yrs. 
In this case, the approach of Ref. [7j performs slightly worse of the exponent expansion for 
n — 1, and slightly better for n — 2. However, also in this case the former breaks down for 
small values of the volatility, generating unphysical transition densities. 
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io- 7 



1 2 
Order of the Expansion (n) 



Figure 3: Accuracy of the exponent expansion for the transition density of the Cox-Ingersoll- 
Ross model (O, for a = 0.0721, b = 0.219, a = 0.06665, and x a = 0.06. Comparison 
between the exponent expansion (stars) and the approximation of Ref. [7] (squares) for 
At = 0.5. The uniform error of the Euler approximation is also reported for comparison. 
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Figure 4: Accuracy of the exponent expansion for the transition density of the Cox-Ingersoll- 
Ross model ([Ml), for a = 0.0721, b = 0.219, a = 0.06665, and x Q = 0.06. Probability density 
function for At = 1.5 for different values of n: n — (dotted), n — 1 (long dashed), n = 2 
(short dashed), n = 3 (continuous), Euler (dot-long dashed), exact (crosses). The inset is 
an enlargement of the region of the maximum. On this scale, the estimates for n = 2 and 
n = 3 still appear coincident. 
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2.3.3 Constant Elasticity of Variance diffusion 

As a last example we consider the Constant Elasticity of Variance model: 

dY t = a(b - Y t )dt + ayYfdWt 



(66) 



Here we consider for brevity only the case p > 1. The transformation to a process with 
constant (unit) variance is X t = 7(Y" t ) = Y^ p /a y (p — 1) and gives, according to Eq. © 



, n Cl _JL_ 

MsW = ^ c 2 x + c^j- 1 

X 



(67) 



with ci = p/2(p - 1), c 2 = a(p - 1), and c 3 = -ab(p - l)P/^-^al /{p 1} . In this case the 
effective potential reads 



Vcff(x) = l;H x (x) 2 



Cl c 2 

2^ + Y + 2(p-l) 



^—Cg^/CP-l) 



(68) 



and the first three terms of the expansion are: 



W (a5, x ) 

M/i(x,x ) 
W 2 (x,a;o) 



ci log ; 



2/ 2(2p-l) 
+ 26(p-l)^a^( 

F(a;) - F(a:o) 



(2p-l)(ar a -^) 



2p-l 

a; p- 1 — at. 



2p-l 
p-1 





1 



2Ax 
1 

2Ax 2 



[y e ff(x) + V cS (x ) - 2W l (x,x Q )] 



(69) 
(70) 
(71) 



ith 



F(z) 



z 3 3p — 1 



c§(p-l) st 



2cic 3 (p-l) _p 2c 2 c 3 (p-l) 3p^2 

2cic 2 zH z»-H — z p- 1 +Vx{z) 

p 3p — 2 



(72) 



Similarly to the examples considered previously, also for the Constant Elasticity of Vari- 
ance model we find a very fast convergence of the exponent expansion for At ~ 1, and 
a performance generally similar to the one of the approach of Ref. [7j, for values of the 
volatility large enough. 



3 PRICING KERNELS OF CONTINGENT CLAIMS 

3.1 Path integral formulation 

The exponent expansion can be generalized to obtain an approximation of the pricing kernels 
of 'standard' derivatives. This can be done by formulating the pricing problem within 
Feynman's path integral framework JSKTU]. Here we indicate as 'standard' any contingent 
claim written on the underlying, Y t , whose value at time t = 0, Vq, can be expressed as an 
expectation value of a functional of a certain type, namely 

V (At,y ) = E\p(Y At )F[Y u ] y ] , (73) 
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where 



F[Y U 



exp 



A/ 



duV F [Y u 



(74) 



and P(Ya*) is a payout function. Above At is the time to expiry, and the expectation 
value is performed with respect to the probability measure defined by the diffusion for Y t 
that we assume of the form ([!}. European Vanilla options, zero coupon bonds, caplets, 
and floorets clearly belong to this family of contingent claims. In addition, other path- 
dependent derivatives, like barrier or Asian options can be expressed in this form (see e.g., 
Refs. Hi). 

Similarly to the case of the transition probability, it is in general convenient to introduce 
an auxiliary diffusion with constant volatility of the form ([5]) by means of the integral 
transformations ([4]). Then, the expectation value in (|T3[) can be transformed in an integral 
over the distribution generated by such auxiliary diffusion by means of the usual Jacobian 
transformation (|10|> . As a result, the value of the option can be in general written as: 



V (At,x ) = E P(X At )F[X u ] 





-L 







dxP(x)K(x,At\x ) 



(75) 



where D is the domain of the auxiliary process as defined by the relative stochastic differ- 
ential equation ([S]), and K(x,At\xa) is the pricing kernel. The latter can be expressed in 
terms of a path integral as it follows [9l [TO] 



with 



®{x,At\ Xo ) 



K(x,At\x ) = < 

c(At)=x 



$(x,At\x Q ) 



At 



du L ff[x(u),x(u)] 



(76) 
(77) 



T>[x(u)] exp 

lx(0)=x o 

where Wo is given by Eq. (|16p and the functional L c ff [x{u) 1 x(u)] is the effective Euclidean 
Lagrangian 

L oB [x(u),x(u)} = —^x(u) 2 + V oS (x) (78) 



2a 



(x(u) = dx(u)/du) with the effective potential, V c f[(x), defined as: 

V c s(x) = — ^^(a;) 2 + \d x {i x (x) + V F (x) 
lot I 



(79) 



Finally, the measure T>[x(u)) is defined by discretizing each path x(u) connecting x(0) = xo 
and x(T) — x. This can be done by dividing the time interval [0, T] into P intervals so that 
x n = x(u n ) (u n = nT/P with n = 0, P), and by integrating the internal P — 1 variables 
x n over the domain D. The path integral J T>[x(u)] is then obtained as the limit for P — > oo 
of this procedure, namely 



f V[x{u)] = lim {2^alAt)- p l 2 ]J [ dx n . 
J p^°° n=1 J D 



(80) 



It is well known form the physical sciences that the path integral $(a;, At|ro) satisfies 
the partial differential equation Eq. (|32|) [12j [13] . Note that this is consistent with the fact 
that, by definition, for Vf(x) = the pricing kernel coincides with the transition density of 
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the underlying diffusion process for X t . In particular, as observed in Sec. 12.2.1] one can use 
Eq. ([32)) to derive the exponent expansion for Q(x, Ax\ro) using the trial form (|34p . As a 
result, the same expressions Eqs. (l23M2"5"j) derived for the transition density hold true also for 
the pricing kernel, provided that the effective potential (|79p replaces the one in Eq. (|26p. 

3.2 Correspondence with Quantum Mechanics 

It is interesting to note that the path integral formulation of the pricing kernel (j77p is 
mathematically equivalent to the quantum mechanical description of the thermodynamic 
properties of an ideal gas of particles moving in the potential hV e s(x) (h is the reduced 
Planck's constant giving the correct energy dimensions). The complete correspondence is 
obtained by identifying a 2 x — ► H/m and At — * H/ksT where m is the mass of the particle, T 
is the temperature, and ks is the Boltzmann constant. With this prescription, $(aj, At\xo) 
becomes the so-called single particle density matrix, and the results of Makri and Miller 
[5] can be readily recovered. In addition, it is straightforward to show using Eqs. ([3171) that 
the exponent expansion of its diagonal elements, &(x Q , At\xg), are consistent with the so 
called Wigner expansion for the high-temperature limit. Finally, performing the analytic 
continuation known as Wick rotation At — > iht allows one to obtain the single-particle 
quantum propagator. In this case (|32p becomes the celebrated Schrodinger equation. 

This correspondence provides an alternative justification of the choice of the exponential 
ansatz in Eq. (|12p . Indeed, this is the form in which can be expressed in general the 
quantum mechanical propagator or the single particle density matrix |12[I13|. Furthermore, 
it has been shown for the quantum problem |22j that the exponent expansion up to third 
order in At and first order in H/m can be derived starting from the short time semiclassical 
propagator obtained through a saddle point analysis of the limit H/m — > [25]. Indeed, it 
can be shown that the second order correction W2(x, xq), Eq. ([24|) . is the so-called van-Vleck 
determinant of the saddle point expansion. This explains why the accuracy of the present 
scheme is preserved in the corresponding regime of low volatility, as it was anticipated in 
Sec l2.2l and illustrated in Sec l2.3l 

3.3 An example: Zero Coupon Bond 

In this Section we illustrate the prescriptions outlined above by applying the exponent 
expansion to the calculation of a zero coupon bond within the Vasicek and Cox-Ingersoll- 
Ross models. The zero coupon bond is a financial instrument that provides at time t = At 
a payout of one unit of a certain notional. Its value at time t = can be expressed therefore 



which is of the standard form given by Eqs. ([73")) and ([7l| . with Vp[X u ] = X u , and P[X& t ] = 
1. 

As a result, the exponent expansion for the kernel Eq. (j76p can be easily derived giving 
for the Vasicek model 



as 




(81) 



W 1 {x,x ) 
W 2 (x,x ) 



W§(x,x ) 
W°(x,x )- 



X + Xq 



(82) 
(83) 



2 



W 3 (x,x ) 



(x-b) 



(84) 
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Figure 5: Accuracy of the exponent expansion for the pricing kernel Eq. (f?6")) of the Vasicek 
model @2J), for a = 0.0717, b = 0.261, and x = 0.1. Symbols as in Fig. gj The panel on 
the left shows the maximum absolute error as a function of the order of the expansion n. 



where W^(x,xq) are the expressions obtained for the transition probability Eqs. (f47|) - (f50|) 
of Sec. 12.3.11 For the Cox-Ingersoll-Ross model instead we get: 



Wi(x,x ) = W?(x,x ) + — (x 2 +x 2 Q + xx ) 



W 2 (x,x ) 



W$(x,x ) 



24 



(85) 
(86) 



with W®(x, Xq) given by Eqs. (fHTTj) . and W^(x, Xq) related to the previous quantities as 
in Eq. flBTJ with 

„1 



(87) 



G°(z) as in Eq. (|62|) . and d z V c f[(z) = d z V® s (z) + <r 2 z/2. The exponent expansion for the 
pricing kernel can be compared with the exact results that can be shown to read for the 
Vasicek model ((42]) 

exp [(a; - x )/a - At(b - a 2 /2a 2 )] 



K ex (x,At\,x ) 



2U/2 



exp 



(27TCT 

[(x -b + a 2 /a 2 ) e 



aAt 



(y-b + a 2 /a 2 ^ 



2a 2 



20 



0123456789 10 

T 



Figure 6: Zero Coupon Bond for the Cox-Ingersoll-Ross model. Parameters and symbols as 
in Figgl 



with a given by Eq. (|4"4")h and 



K ex (x,At\,x ) = -exp 



-^x 2 -x 2 ) + (2ab-^)logf- 



^a 2 bAt/a 2 



2a 2 sinh [jAt/2] 



exp 



4a 2 



(x 2 + xl)coth[jAt/2] 



xx "/ 



q V2cr 2 sinh[7A</2] 



(89) 



with 7 = V» 2 + 2cr 2 , for the Cox-Ingersoll-Ross one. As illustrated in FigJSl similarly 
to the case of the transition probability, the exponent expansion provides a remarkably 
good, and fast converging approximation of the exact pricing kernel for financially sensible 
parameterizations, and for a sizable value of the time step At. 

Finally, the zero coupon bond can be obtained by numerical integration of the pricing 
kernel according to Eq. (|75p . The corresponding results are shown in Fig. [6] confirming 
once more the quality of the approximation. In a similar fashion, one can obtain systematic 
approximations for caplets, floorets, and other simple interest rate derivatives whose value 
depends only on the value of the instantaneous short rate at time t = At. It is also possible 
to generalize this approach to path dependent contingent claims, like Asian options. 



4 EXTENDING THE TIME STEP: PATH INTEGRAL 
MONTE CARLO METHODS 

For an extended time interval T, the calculation of the transition density or, more in general, 
of the pricing kernel (|76|) can be performed by discretizing the path integral (|77| , and taking 
the limit of large number of time slices (P — > oo) according to the standard Trotter product 
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formula: 



K(x,At\x ) ~ (2TT<rlAt) 



-P/2 e -W a {x,x ) 



P-1 

\\ J dx n 



x exp 



k=P 



2alAt 



k=l 



71=1 

k=P 



k=l 



(90) 



with At = T/P, ip = x, and the effective potential given by Eq. (|79j) . It is worth re- 
marking that, for the case of the transition density, by interpreting the latter equation 
as the Chapman-Kolmogorov property of Markov processes [TT] one obtains the following 
approximation of the short-time propagator 



^Trottor(£, At\x Q ) 



g— W {x,x a ) 



exp 



(x - Xp) 2 
2alAt 



^-(V c b(x) + V e s(xo] 



(91) 



However, in contrast to the n = 1 exponent expansion, the latter expression is not strictly 
correct up to order At, and only in the limit P — > oo the difference becomes negligible. 

In general, to obtain an accurate result for the pricing kernel for an extended time period 
T one has to increase the number of time slices, or Trotter number P, until convergence is 
achieved. By replacing the Trotter formula with the improved short-time kernel obtained 
through the exponent expansion (|12[) one achieves a faster convergence with the Trotter 
number, thus significantly reducing the computational burden. In this case the finite-time 
expression of the pricing kernel reads 

p-i 



K(x,T\x ) ~ {2na 2 x At)- p / 2 JJ f dx n 

n=l •' D 



x exp 



2 At. ^ V 



2alAt 



fe=i 



fc=i 



(92) 



with W(x, xo, At) given by Eq. ([T5)l . 

Equation (|92[) allows one to obtain the transition density or the pricing kernel for a 
standard derivative through the calculation of a multidimensional integral over the variables 
Xi, . . . , Xp—i. The latter integration is ideally suited for Monte Carlo methods either in real, 
or in Fourier space [26] . the specific choice depending on the particular problem at hand. 
In addition, importance-sampling schemes, e.g., by means of the Metropolis algorithm [27j . 
can be easily applied in order to reduce the computation time. However, in order not to 
introduce a systematic bias in the result a particular attention has to be paid in order to 
sample accurately the configuration space. 

The most straightforward way to perform a Monte Carlo quadrature of Eq. (|92|) , is to 
realize that a simple Markov chain x = (x±, . . . , Xp—\) 



x„__i + <7 x VAtZ n , 



(93) 



with Z n , sampled from a standard normal distribution, generates an ensemble of walkers 
distributed according to 



p(x 1 ,...,x P - 1 \x ) = (2ira 2 x At)-( p - 1 » 2 



exp 



1 P_1 

— ^{Xk-Xk^f 



2alAt 



fc=i 



(94) 
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As a result, the pricing kernel (f9"2"|) can be obtained as the average over the random paths 
generated according to Eq. ([93]) of the following estimator: 



C(x, x P = x) = {2na 2 x At)- 1/2 exp 



l p 
j— (x- zp-i) 2 -22,W{x k ,x k -i) 



2alAt 



k=i 



(95) 



A remarkable property of the path integral approach is that K(x, At\xo) for any final point 
x can be evaluated with a single Monte Carlo simulation by appropriately reweighting the 
accumulated estimator. In fact the distribution of walkers p(xi, . . . , xp-i \xq) is independent 
of the final point so that K(x', At\x ) can be calculated by averaging 0(x, xp — x'). In 
addition, the latter quantity can be efficiently obtained by means of the following reweighting 
procedure 

0(x,x P = x') = 0{x,x P = x) ^f ,X } (96) 



with 



W(x,x) = exp 



1 



2alAt 



W(x,x) 



{x- xp^if -W^tXp-x) 



(97) 



Expectation values of the form (|75|) on a time horizon T can be calculated by integrating 
over the final variable giving: 

p 

V {T,x ) = [ dxP(x)K{x,T\x ) ~ (2^At)- p / 2 "TT / dx n P{x P ) 
Jd ; Jd 



exp 



1 P P 

2-TT ^2{%k - Xk-l) 2 - W(x k ,X k -l) 



2alAt 



fe=l 



(98) 



This can be simulated by extending the Markov chain ((93)) up to step x — xp, and accumu- 
lating the estimator 



■p(x, xp = x) = P(x) exp 



,X k -l) 



fc=l 



(99) 



Remarkably, within the path integral approach, the sensitivities of such expectation values 
with respect to the model parameters (the so-called Greeks) can be computed in the same 
Monte Carlo simulation, thus avoiding any numerical differentiation. Indeed, indicating 
with 9 a generic parameter, under quite general conditions [28], one has 



d e Vo(T,x ,e) = / dx[K e (x,T\xo)deP(x,0)+P(x,0)d e K e (x,T\xo)] ■ (100) 

JD 

As a result the sensitivity dgVo(T, Xq, 9) can be calculated by means of the estimator: 



Q(x, xp = x) = exp 



,Xk-l) 



k=l 



{dgP + Pde log Kg) . 



(101) 



Higher order sensitivities can be obtained in a similar fashion. 

The convergence with the Trotter number P of the path integral Monte Carlo estimates is 
illustrated in Fig. [7] for the calculation of the first five moments of the T = 40 yrs transition 
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Figure 7: Convergence with the Trotter number P of the normalization, and of the first 4 
moments of the 40 years transition density for the Cox-Ingerssol-Ross model. Parameters as 
in Fig[3] In the first top panels the relative error with respect to the exact result is reported. 
Lines are quadratic fits. 



probability of the Cox-Ingersoll-Ross model (p>4"l) . The finite P estimates converge very 
rapidly with 1/P. In particular, for the case considered, P = 20 already provides estimates in 
agreement with the exact result within statistical uncertainties. In general, as also shown in 
the figure, a convenient indicator of the convergence is the zeroth-moment or normalization 
of the distribution. The calculation of this quantity allows in general to assess the level 
of convergence without performing a complete scaling with P, thus saving computational 
time. 

5 EXTENDING THE TIME STEP: GREEN'S FUNC- 
TION BACKWARD INDUCTION 

The Green's function backward induction is a deterministic approach alternative to partial 
differential equation methods, introduced in financial context by Rosa-Clot and Taddei [ID] . 
This is based on the path integral representation of standard contingent claims (|73j) . that 
we rewrite here for convenience as 

E[G T {Y U ) I Y ] . (102) 
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Here Gt is the payout functional dependent on the realization of the path of the underlying 
up to expiry T, and Y t is the general process following the diffusion fl}. As partial differential 
methods, the Green's function backward induction is efficient for low-dimensional problems 
but allows one to include naturally early exercise features in the contingent claim, i.e., to 
treat expectation values of the form 

maxE[G T {Y u )} (103) 

T 

where < r < T is the stopping time of the process [TT] . 

For standard contingent claims, Eqs. (|73[) and (f?4")) . the functional Gt[Xi] can be dis- 
cretized in the form, 

JV 

G T [Y u ]^l[g^(Y t ). (104) 

i=0 

and the conditional expectation value (|102|) is given by 

„ „ N N 
E[G T [Y U ] I Y ] ~ /... / l[dYi g( N \Y N ) \{p{Y 3 ^ \ ^-i.tj-i), (105) 

i=l j=X 

where the function, p, is 

~p{Y u ti | Y i _ 1 ,t i _ 1 ) = pyiY^U - u-i I g^HYi-x) , (106) 

and p y (Y i} ti — tj_i \ Yi-x) is the transition probability associated with the process ([J). This 
is the Green's function of the partial differential equation associated with the calculation of 
the expectation value (|102[) , hence the name of the method. 
Let us now consider a single integration 

J<lY l W^\Y.,t l )W,U\Y.- 1 ,t.- 1 ). (107) 

If we approximate this integral by using a numerical quadrature rule, we obtain the following 
algebraic relation 

M 
7=1 

where the matrices, p W , are dehned by 

Pafi = P( z a,U+l | Zp,U), (109) 

the quantities, w a , and z a , are the weights and the grid points, respectively, associated with 
the integration rule, and a, (3, 7 = 1,...,M. In conclusion, the expression (|105|1 can be 
written as 

M 

E[G T [Y(r)] I Y = z a ] * £ G« 1 . . . G^Z^AZ^l tif . (^0) 

7i>— i7w=l 

where G^p = w aPaP' an< ^ ffw = 9^ N \ z in)- Therefore, we have reduced the evaluation 
of the expectation value of a functional to the product of ./V matrices with dimension M. 
By starting the calculation from the right (backward induction), we need to memorize just 
linear arrays, while the matrix elements, G^, can be computed step by step. In practice, 
one can follow the following algorithm: 



25 



i. Set u a = g a ' (a = 1, . . . , M), and i = N - 1. 



M 



ii. Set v a = 



(«=l,...,M). 



/3=1 



iii. If z > then set u a — v a (a = 1, . . . , M), i = i — 1, and go to ii. 

Here the arrays, u Q , and are two working vectors. 

The inclusion of early exercise features is completely straightforward and simply involves 
the following steps: 



(a=l,...,M). 

iii. If i > then set u a — v a (a = 1, . . . , M), i = i — 1, and go to ii. 

The only difference is in the point ii, where we have now a test operation, and the function 
f(z a ,ti) is the exercise value of the option. 

The exponent expansion can significatively speed up the Green's function backward 
induction presented above. In fact, as illustrated in the previous Sections, it allows one 
obtain reliable approximations of p(Yi,ti | Yi_i,ti_i) for a time step (U — much 
larger than those that can be safely used with a simple Euler discretization. As a result, 
one can reduce the discretization bias to acceptable levels with a much smaller number of 
intermediate time steps AT. The application of this procedure to local volatility models is 
currently in progress and will be presented elsewhere 29J. 



Closed-form approximations of non-liner diffusions are of primary importance in a variety 
of fields of quantitative Finance. In Econometrics, they are crucial for an efficient maximum 
likelihood estimation of the parameters of continuous time processes. In derivative pricing, 
they allow one to develop effective approximation schemes or to improve the efficiency of 
numerical approaches. 

In this paper we have presented an effective method to produce a family of closed-form 
approximations of the transition probability of a general diffusion process. Such approxima- 
tion, dubbed exponent expansion, is based on an exponential ansatz of the transition prob- 
ability for a finite time interval At, and a series expansion of the deviation of its logarithm 
from that of a Gaussian distribution. Through this procedure the transition probability 
is obtained as a power series in At which becomes asymptotically exact if an increasing 
number of terms is included, and provides remarkably accurate results even when truncated 
to the first few (say 3) terms. This approach can be easily implemented, and involves the 
calculation of simple one dimensional integrals. In addition, it applies to a very wide class 
of volatility functions, even to those that do not have a simple analytic expression but they 
are specified through a numerical interpolation. 



i. Set u a — g a (a — 1, ... , M), and i = N — 1. 

ii. Set mi'' = f(z a , ti), and v a — max I ~S~^ up G' 




f3a 1 



W, 




6 CONCLUSIONS 
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We have shown that the exponent expansion produces very accurate results for inte- 
grable diffusions of financial interest, like the Vasicek and the Cox-Ingersoll-Ross models. 
In particular, we have compared our results with those obtained with the state-of-the-art 
approximation of discretely sampled diffusions [7j, that shares a similar rationale. We find 
that the exponent expansion provides a similar accuracy in most of the cases but it is 
more stable in the low-volatility regimes. Furthermore the implementation of the exponent 
expansion turns out to be simpler. 

By introducing a path integral framework we have generalized the exponent expansion 
to the calculation of the pricing kernels of financial derivatives, and we have shown how 
to obtain accurate approximations for the price of simple contingent claims. We have also 
shown how the exponent expansion can be naturally used to increase the efficiency of both 
deterministic and stochastic numerical simulations. A systematic study of the efficiency of 
this approach for the pricing of exotic derivatives, and the calibration of local and stochastic 
volatility models will be the object of future investigations. 
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